Designing peptides predicted to bind to the omicron variant better than ACE2 via computational protein design and molecular dynamics

Brought about by severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2), coronavirus disease (COVID-19) pandemic has resulted in large numbers of worldwide deaths and cases. Several SARS-CoV-2 variants have evolved, and Omicron (B.1.1.529) was one of the important variants of concern. It gets inside human cells by using its S1 subunit’s receptor-binding domain (SARS-CoV-2-RBD) to bind to Angiotensin-converting enzyme 2 receptor’s peptidase domain (ACE2-PD). Using peptides to inhibit binding interactions (BIs) between ACE2-PD and SARS-CoV-2-RBD is one of promising COVID-19 therapies. Employing computational protein design (CPD) as well as molecular dynamics (MD), this study used ACE2-PD’s α1 helix to generate novel 25-mer peptide binders (SPB25) of Omicron RBD that have predicted binding affinities (ΔGbind (MM‑GBSA)) better than ACE2 by increasing favorable BIs between SPB25 and the conserved residues of RBD. Results from MD and the MM-GBSA method identified two best designed peptides (SPB25T7L/K11A and SPB25T7L/K11L with ΔGbind (MM‑GBSA) of −92.4 ± 0.4 and −95.7 ± 0.5 kcal/mol, respectively) that have better ΔGbind (MM‑GBSA) to Omicron RBD than ACE2 (−87.9 ± 0.5 kcal/mol) and SPB25 (−71.6 ± 0.5 kcal/mol). Additionally, they were predicted to have slightly higher stabilities, based on their percent helicities in water, than SBP1 (the experimentally proven inhibitor of SARS-CoV-2-RBD). Our two best designed SPB25s are promising candidates as omicron variant inhibitors.

To inhibit the virus from entering human cells, blocking binding interactions (BIs) between ACE2-PD and RBD using peptide inhibitors, neutralizing antibodies and small-molecule drugs have been extensively investigated [10][11][12][13][14][15][16].Due to their high structural compatibility with the surface of a protein target, and their abilities to disrupt protein-protein interaction, peptides can be utilized as inhibitors that disrupt ACE2-PD and RBD binding [17,18].
Computational approaches including molecular dynamics (MD) and computational protein design (CPD) have been utilized for designing peptide binders of RBD of the wuhan variant and other variants [19][20][21][22].To inhibit ACE2 and RBD binding, CPD and MD were employed in our previous studies to design novel 25-mer peptide binders (SPB25), which were derived from ACE2-PD's α1 helix, that have better predicted SARS-CoV-2-RBD binding affinities (ΔG bind (MM-GBSA) ) than ACE2 and 23-mer peptide binder (SBP1) [23,24], which is the experimentally proven inhibitor of SARS-CoV-2-RBD [25].Nonetheless, the knowledge on BIs between SPB25 and Omicron RBD as well as on whether our strategy, increasing favorable BIs between SPB25 and the conserved residues of RBD, can be used to create novel 25-mer peptides that have predicted ΔG bind (MM-GBSA) to Omicron RBD better than ACE2 is limited.
Based on ACE2-PD's α1 helix (residues , this study used CPD (Rosetta) as well as MD (AMBER) to generate novel SPB25 that have better predicted ΔG bind (MM-GBSA) to Omicron RBD than ACE2.We utilized the design strategy of increasing favorable BIs between SPB25 and the conserved residues of RBD.The designed SPB25 that have predicted ΔG bind (MM-GBSA) to Omicron RBD better than ACE2 are promising peptides that may be utilized to inhibit BI between ACE2 and Omicron RBD.

Computational design of novel SPB25s
In this study, the design template for CPD (Rosetta [29]) was the structure of SPB25 in complex with Omicron RBD.This work employed the design strategy of increasing favorable BIs between the conserved residues of RBD (Y421, L455, F456, G485, F486 and Y489) [30] and residues 21-45 of ACE2 so that the designed SPB25 have better ΔG bind (MM-GBSA) to Omicron RBD than ACE2.Designed positions were selected if favorable BIs could potentially be formed upon mutations between RBD's conserved residues and the side chains of designed positions.Employing CoupledMoves protocol [31,32] of RosettaDesign module (Rosetta3.11)and beta_nov16 energy function, structures of selected designed residues/positions were redesigned, repacked and minimized.Standard amino acids except G and P were allowed in each designed position.Additionally, residues within 10 Å of each designed position were repacked and energetically minimized.400 runs were independently conducted for each design, producing 400 conformations of designed sequences, where some sequences may have various conformations.To calculate ΔG bind (Rosetta) in Rosetta Energy Unit (REU) of each conformation, an interface analyzer [33,34] module (Rosetta3.11)was used.ΔΔG bind (Rosetta) upon mutation was determined by deducting ΔG bind (Rosetta) of SPB25 from that of the designed conformation/sequence.MD was performed on designed conformations of all designed positions with ΔΔG bind (Rosetta) < 0 REU, and the molecular mechanics-generalized born surface area (MM-GBSA) method [35][36][37] was utilized to determine if their ΔG bind (MM-GBSA) values are better than that of SPB25.The single mutants with better predicted ΔG bind (MM-GBSA) than SPB25 were used to generate the double mutants of SPB25, and subsequently simulated to determine if their ΔG bind (MM-GBSA) values are superior to that of SPB25.

MD and analyses
Complex structures of Omicron RBD binding to ACE2 and designed peptides were built in isomeric truncated octahedral boxes of TIP3P water utilizing the buffer distance of 13 Å as well as force field parameters of protein.ff14SB[38] and GLYCAM06j-1 [39] in AMBER18 [28].To remove interactions that are unfavorable, all systems were energetically minimized, using the five-step procedure [23,24,40].Employing different restraints on proteins, each minimization step contains 2,500 steps of steepest descent and 2,500 steps of conjugate gradient.In the first step, hydrogen atoms and water molecules were energetically minimized, while restraining proteins' heavy atoms using a force constant of 10 kcal/(mol Å 2 ).In the second, third and fourth steps of minimizations, force constants of 10, 5 and 1 kcal/(mol Å 2 ) were subsequently applied to restrain proteins' backbones, respectively.Lastly, the whole system was minimized without any restraining force.Subsequently, the GPU (CUDA) version of PMEMD module [41][42][43] was used for MD with the periodic boundary condition.To constrain all bonds with hydrogen atoms to allow the time step of 0.002 ps, SHAKE [44] was employed.To control each system's temperature, the Langevin dynamics method with a collision frequency of 1.0 ps -1 was used.For heating, the temperature of each system was increased from 0 K to 310 K in the NVT ensemble during 200 ps MD, while the backbones of the proteins were restrained with the force constant of 10 kcal/(mol Å 2 ).Each system was additionally equilibrated in the NVT ensemble without any restraining force at 310 K for 300 ps.In the production run, each system was subject to MD in the NPT ensemble at 310 K and 1 atm for 100 ns.
The structural stabilities of all systems were measured based on their Root Mean Square Deviation (RMSD) values with respect to their minimized structures.Further analyses were conducted on 80-100 ns trajectories, which have stable RMSD values.To predict each complex's binding affinity, the MM-GBSA technique was utilized to compute their total binding free energies [ΔG bind (MM-GBSA) ].The MM-GBSA method uses the following equation to compute ΔG bind (MM-GBSA) .
where ΔG bind is the total binding free energy of the system that is defined as the difference between the free energies of the complex (G complex ), the receptor (G receptor ), and the ligand (G ligand ).ΔG bind is calculated as the sum of the molecular mechanics free energy (ΔE MM ), the solvation free energy consisting of polar (ΔG polar(GB) ) and nonpolar contributions (ΔG non-polar (SA) ), and the entropy (-TΔS) in the gas phase.The ΔE MM term comprises ΔE vdw (van der Waals) and ΔE ele (electrostatic).The polar and non-polar contributions are estimated by the Generalized-Born (GB) implicit solvation model and the molecular solvent accessible surface area (SASA).In this work, the entropic contribution (-TΔS) was not included in the calculation of ΔG bind because the nmode module of AMBER predicts this term with high computational cost but not high accuracy [45,46].Designed SPB25s with ΔG bind (MM-GBSA) better than ACE2 were selected for additional analyses including BIs and per-residue free energy decomposition (PFED).Hydrogen-bond (H-bond) occupations were determined to identify H-bonds of all systems.The following rules were employed to consider an occurrence of a H-bond: (i) a proton donor-acceptor distance �3.5 Å and (ii) a donor-H-acceptor bond angle �120˚.In this work, three levels of occupations were defined: (i) strong H-bonds with H-bond occupations > 75%, (ii) medium Hbonds with 75% � H-bond occupations > 50%, and (iii) weak H-bond interactions with 50% � H-bond occupations > 25% [23,24,40,47].To measure the peptide helicity of each system, Define Secondary Structure of Protein (DSSP) was utilized to compute the percent helicity from the summation of the percentage of α-, π-and 3-10-helix structures [48].

MD and calculations of binding free energies
Complex structures of all 16 designed 25-mer peptides, SPB25 and ACE2 binding to Omicron RBD were simulated for 100 ns.Using the MM-GBSA technique, ΔG bind (MM-GBSA) of all complexes were computed to assess whether ΔG bind (MM-GBSA) of designed SPB25s are better than that of SPB25.To determine the stabilities of each system, Root Mean Square Deviation (RMSD) of all and backbone atoms were computed (S1 Fig) .RMSD plots show that it likely took around 80 ns for each system to reach equilibrium.As a result, further analyses were conducted on the 80-100 ns trajectory of each system.
The binding free energy components of designed SPB25s with better ΔG bind (MM-GBSA) than ACE2 The binding free energy components of ACE2, SPB25 and two designed SPB25s with better predicted ΔG bind (MM-GBSA) to Omicron RBD than ACE2 are shown in Fig 2 .The main contributor to the favorable ΔG bind (MM-GBSA) to Omicron RBD of ACE2, SPB25, SPB25 T7L/K11A , and SPB25 T7L/K11L is the electrostatic interaction term.Other terms that also favorably contribute to ΔG bind (MM-GBSA) are the non-polar solvation and van der Waals energy terms.The term that unfavorably contributes to ΔG bind (MM-GBSA) is the polar solvation term.SPB25 T7L/K11A and SPB25 T7L/K11L are the designed 25-mer peptides that have the best ΔG bind (MM-GBSA) to Omicron RBD (ΔG bind (MM-GBSA) values of −92.4 ± 0.4 and −95.7 ± 0.5 kcal/mol, respectively).Their ΔG bind (MM-GBSA) are better than that of SPB25 by 20.8 ± 0.6 and 24.1 ± 0.7 kcal/mol, respectively, and better than that of ACE2 by 4.5 ± 0.6 and 7.8 ± 0.7 kcal/ mol, respectively.As compared to those of SPB25, the major contributions to the favorable BI of SPB25 T7L/K11A and SPB25 T7L/K11L to Omicron RBD are the increases in their favorable electrostatic interaction, van der Waals energy, and non-polar solvation terms.Nevertheless, their unfavorable polar solvation terms are higher than that of SPB25.
Fig 4 illustrates key predicted BIs between the two best designed 25-mer peptides and Omicron RBD.The overall binding poses of these two designed SPB25s and ACE2 binding to omicron RBD are relatively similar.In terms of predicted interactions to Omicron RBD, the total no. of H-bonds of SPB25 T7L/K11A is higher than those of ACE2 and SPB25.Additionally, no. of strong H-bonds of SPB25 T7L/K11A is higher than those of ACE2 and SPB25, supporting the result that its ΔG bind (MM-GBSA) is better than ACE2 and SPB25.D18 of SPB25 T7L/K11A had four strong H-bonds with R498 and Y501 of Omicron RBD.Additionally, H14 of SPB25 T7L/K11A had three strong H-bonds with R493 and Y501 of Omicron RBD.D10 and E15 of SPB25 T7L/ K11A had four medium H-bonds with R403 and R493 of Omicron RBD.Furthermore, other residues such as Q4 and Y21 had H-bonds with Omicron RBD.In terms of π interactions between SPB25 T7L/K11A and Omicron RBD, no. of π interactions of SPB25 T7L/K11A is more than that of SPB25.H14 of SPB25 T7L/K11A had two cation-π interactions with R403 and R493 of Omicron RBD.In addition, there were five π-π interactions (F8 -Y489, H14 -Y501, H14 - H505, Y21 -Y501 and Y21 -H505), one σ-π interaction (Y21 -Y501@HA), and two alkyl-π interactions (A5 -Y489 and H14 -R493) between SPB25 T7L/K11A and Omicron RBD.
https://doi.org/10.1371/journal.pone.0292589.t002Peptide helicities of designed SPB25s with better ΔG bind (MM-GBSA) than ACE2 Figs 5 and S2 show the percent helicities and structural stabilities (RMSD plots) in water of designed 25-mer peptides with better predicted bind affinities to Omicron RBD than ACE2, respectively.The percent helicities of the N-and C-termini of the two best designed peptides are lower than those of middle residues due to their high flexibilities.Overall trends of percent helicities of SPB25 T7L/K11A and SPB25 T7L/K11L are slightly higher than those of SBP1.

Discussion
The omicron variant (B.1.1.529)was the important variant of concern that was responsible for the COVID-19 pandemic.Similar to other variants, Omicron RBD firstly attaches to ACE2-PD to enter human cells.Disrupting BIs between ACE2-PD and Omicron RBD to prevent coronavirus from infecting and destroying human cells is a promising COVID-19 therapy.Since peptides have more similar interactions to native protein-protein interactions and functional groups than small molecules, which can be ineffective in disrupting large proteinbinding interfaces [17,18], peptides can be employed as inhibitors of SARS-CoV-2 to inhibit protein-protein interactions at their binding interfaces.
To design novel 25-mer peptides with high potential to bind to Omicron RBD better than ACE2, we employed CPD, using the residues 21-45 of ACE2-PD's α1 helix as a template, and MD.The design strategy of this study was increasing favorable BIs between SPB25 and the conserved residues of RBD (Y421, L455, F456, G485, F486 and Y489).Q4 (24), T7 (27), F8 (28), D10 (30), K11 (31) and H14 (34) were chosen as designed positions because their side chains are in the orientations that can possibly form favorable BIs with Omicron RBD upon mutations.Standard amino acids, except G and P, were allowed for all designed positions.After CPD by Rosetta, 52 designed SPB25s that have single mutations were generated.Sixteen designed SPB25s with superior ΔG bind (Rosetta) to SPB25 (ΔΔG bind (Rosetta) < 0 REU) were simulated, and their ΔG bind (MM-GBSA) were computed by the MM-GBSA technique and compared to those of SPB25 and ACE2.The predicted binding affinity of ACE2 to Omicron RBD (ΔG bind (MM-GBSA) = −87.9± 0.5 kcal/mol) is better than that of SARS-CoV-2-RBD (ΔG bind (MM-GBSA) = −71.2± 0.4 kcal/mol) [23], supporting the experimental findings that ACE2 bound to Omicron RBD (K D = 38.9 ± 10.5 nM) with higher affinity than the wild type (K D = 75.5 ± 2.1 nM) [9].Additionally, the predicted binding affinity of SPB25 to Omicron RBD (−71.6 ± 0.5 kcal/mol) is better than that of SARS-CoV-2-RBD (−60.3±0.4 kcal/mol) [23], suggesting that it should be able to experimentally bind to Omicron RBD better than that of the wild type.However, the predicted binding affinity to Omicron RBD of SPB25 is worse than that of ACE2 probably because no. of IBRs of ACE2 to Omicron RBD is significantly more than SPB25, which also includes those in the α2-helix and the linker of the β3-and β4-sheets.
The trends of percent helicities in water of SPB25 T7L/K11A and SPB25 T7L/K11L are relatively similar to SPB25 [23] and slightly higher than that of SBP1 [23].Our findings suggest that their stabilities in water may be relatively similar to SPB25 and slightly higher than that of SBP1, and the stabilities of these designed SPB25 should be sufficient for their use as peptide binders.
Using the combination of CPD and MD, we designed promising SPB25s with better ΔG bind (MM-GBSA) to Omicron RBD than human ACE2 receptor and SPB25 by increasing favorable BIs between peptides and the conserved residues of RBD.The results from MD and the MM-GBSA calculation show that the values of ΔG bind (MM-GBSA) of two best designed peptides (SPB25 T7L/K11A and SPB25 T7L/K11L ) are better than those of ACE2 and SPB25.Moreover, their predicted helicities in water are slightly higher than that of SBP1, suggesting that their stabilities are higher than that of SBP1.SPB25 T7L/K11A and SPB25 T7L/K11L are promising peptide candidates that could possibly be utilized to disrupt BIs between Omicron RBD and ACE2.

Fig 1 .
Fig 1.The design template structure of SPB25 binding to Omicron RBD.Omicron RBD and SPB25 are displayed in cyan and orange, respectively.Conserved residues of RBD and the designed positions are labelled in blue and red, respectively.https://doi.org/10.1371/journal.pone.0292589.g001